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ABSTRACT 

We review in brief the development and implementation of the Star integral, a tool 
yielding measurements of correlations much superior to conventional methods. A ver- 
sion for use in pion interferometry is explained. We also show how effects of non- 
poissonian overall multiplicity distributions may be eliminated if desired and quote 
results eliminating statistical biases arising in correlation measurements within small 
samples. 
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1 Point distributions 



Distributions are fundamental to almost any branch of the exact sciences. A point distri- 
bution is characterized by the fact that the object under scrutiny has no intrinsic structure 
or size; it is a point rather than a field. The coordinates of this point may be discrete or 
continuous, real or complex, embedded in a single- or multidimensional space. Physically, 
whenever the objects ( "particles") have a size small in comparison with the embedding space, 
the assumption of a point distribution is justified. 

Typical examples of point distributions are galaxies in the sky and pions in phase space; 
in these cases, the embedding space is continuous. If Xi are the coordinates of iV particles 
measured in a particular "event" (as measured by the detector, a region of the sky, a throw 
of N dice, . . . ), the number of such particles at the point x is 

p x {x) = Y,8{x-X h ). (1) 

ii=l 

In general, the simultaneous behavior ( = correlation) of q of these particles for a given event 
a is given by 

TV 

%(x lt ...,x q )= £ 5( Xl - XI) S(x 2 - XI ) ■ ■ ■ S(x q - X a iq ) . (2) 
=1 

Note that p q is a nonnegative integer while the coordinates X are continuous. Meaningful 
results are extracted by averaging over many events, to yield the g-tuple density 

TV CV 

p 9 = (^)=TV- 1 ^^, (3) 

a=l 

which is pseudocontinuous (actually, a rational number with denominator N ev ). Allowing 
the variables to range over some restricted domain Q, one measures moments of point dis- 
tributions 

fa(fi) = J p q (x 1 ,...,x 9 )dx 1 ...dx 9 = (n [q] ) Q (4) 

which contain information on the correlation strength in this particular region Q. Here 
n' 9 ' = n!/(n — q)\ = n(n—l) ■ ■ ■ (n—q+1), so that £ 9 is actually a factorial moment of the 
distribution. 

Below, we explore one particular type of domain Q which leads to the so-called Star 
integral. The reason for this particular choice of Q is that it maximizes the amount of 
information extracted from a given sample of events while restricting the amount of computer 
time to a minimum JT|. 
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2 Star integral moments 



Conventional measurements of correlations proceed first to discretize the continuous variable 
X ("binning the data") and then to find £ 9 by averaging over all events the counts nj* in 
every bin 0, 




By contrast, the Star integral belongs to the class of correlation integrals, where distances 
between pairs of points X ili2 = \X^ — X i2 \ are computed directly, before binning ||. For 
the Star integral, the domain Q is given by the sum of all spheres of radius e centered around 
each of the N particles in the event. The number of particles ("sphere count") within each 
of these spheres is, not counting the particle at the center X ix , 

N 

h{X h ,e) = ®( e -^M 2 ); «2^«i, (6) 

12 = 1 

and the factorial moment of order q is 

^(e) = (£h(X il ,e)^y (7) 

The obvious similarity to the conventional moment of Eq. (0) should not disguise the fact 
that J2n is a sum over particles rather than bins. 

The Star factorial moment of Eq. (|7|) can be derived rigorously [jy from Eq. @ using for 
Q the equivalent implicit definition 

£f "(e) = / Pg(xi, ...,x q ) 12 6i 3 . . . e lq dx 1 . . . dx q (8) 

with the theta functions 0y = 0(e — \x\ — Xj\) restricting all q—1 coordinates Xj to within 
a distance e of X\. 

Correlation measurements can be made either as a function of the sphere size e — useful 
in looking for self-similarity and fractal structure - - as a function of the distance from 
a fixed center coordinate, as has been the case traditionally. This is implemented in the 
"differentials" of Section 3 below. 

The superiority of the Star integral moments of Eq. ([8]) over the conventional ^ onv arises 
because the artificial discretization inherent in the latter has two bad effects ||. First, it 
leads to instabilities in the measured quantities as the bin size is varied and, second, 
it often sorts particles into different bins when in fact they are quite close together, thus 
unwittingly throwing away information. By contrast, ^ tar is much more stable and has 
smaller errors, especially if the coordinates in use live in a two- or three-dimensional space. 

In the literature on galaxy correlations Q and in the characterization of strange attractors 
0, an approach similar to our Star integral has been used for some time, utilizing, however, 
not the factorial but the ordinary form 

(^pp^r 1 )- (9) 
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We believe that this is an ad hoc approximation to the exact factorial form of Eq. (|7|), 
necessarily breaking down when h becomes small. Current measurements in these fields may 
therefore be suffering from distortion at small e. 

In order to eliminate, among other things, the overall multiplicity, it has become custom- 
ary in high energy physics to measure normalized factorial moments . The denominator 
used for such normalization should be made up of the uncorrelated background, p\. While 
it can be implemented in a number of ways, we prefer the "vertical" normalization, in which 
p\ is integrated over exactly the same domain Q as the inclusive density p q in the numerator. 
Thus for the Star integral, the normalized moment is 



£f ar f p q (x!,..., Xg) 612Q13 ■ ■ ■ Big dxi . . . dx q 
q VV ~ C orm ~ I Pi{x 1 )...p 1 {x q )Q l2 Q 13 ...Q lq dx l ...dx q 



HStar/ \ _ *sg 



We have shown rigorously that the denominator £° orm is given by the following double 
event average: with X°^ i2 = jX? — X i \ measuring the distance between two particles taken 
from different events a and b, 

C rm ( e ) = (e(e^-^5 2 )) 9 ) = (e^W^c)) 9 " 1 ), (ii) 

where the outer event average and sum over %\ are taken over the center particle taken from 
event a, each of which is used as the center of sphere counts h\,{X\. e) taken over other 
events b in the inner event average. We thus see the natural emergence of the heuristic 
procedure of normalization known as "event mixing" [p], [J. 

Apart from the double event average and the appearance of the ordinary power q—1 
rather than the factorial power [q— 1], the similarities between the numerator Eq. (|7D and 
denominator Eq. ( |TT1) are obvious. Both do sphere counts around a given center particle 
at Xi x ] the numerator £^ tar does so within the same event a, while the denominator £° orm 
inserts this center particle into all other events b to perform a similar count there. This is 
shown schematically in Figure 1. 



3 Cumulants and differentials 

Cumulants are combinations of correlation functions constructed in such a way as to become 
zero whenever any one or more of the points x becomes statistically independent of the 
others. This is done so as to strip away the combinatorial background from the correlation 
measurements, 

C 2 (x 1 ,x 2 ) = p 2 (xi,x 2 ) - p 1 (x 1 )pi(x 2 ) , (12) 
C 3 (x h x 2 ,x 3 ) = p 3 (x 1 ,x 2 ,x 3 ) - p^xijp^x^x-i) - Pi(x2)p 2 (x 3 ,x 1 ) 

- Pi{x 3 )p 2 (x 1 ,x 2 )+ 2p 1 (x 1 )p 1 (x 2 )p 1 (x 3 ) etc. (13) 
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event a 



event b 



Figure 1: Sphere counts. On the left is shown a typical event a, with the particles 
denoted as circles. For each particle i\ of a, count all other particles within a sphere 
of radius e; this yields r^-X^e) of Eq. (||) used in the numerator of the Star integral 
moments and cumulants. On the left is shown a different event b, with particles denoted 
as squares. For normalization and cumulants, the same center particle is inserted at 
X j x into event b and a count performed to yield n^X^ , e). Performing an event average 
over all 6-events, one obtains the normalization £n° rm as in Eq. ( |Tl| ) and cumulants as 
in Eq. (pi). 



Integrating these over the Star integral domain, we can find the normalized cumulants 

K^(e) = -^-, (II) 

with 



f q (e) = J C q (x u ...,x q ) ©12613 ...Qigdxx... dx q (15) 

the unnormalized (Star) factorial cumulant. The latter can be written entirely in terms of 
the sphere counts introduced previously; writing in shorthand 

a = £e(e-X-) = n(X?,e), j^i (16) 

j 

b = Ee(e-lJ)=M4e) (17) 

j 

and defining for convenience the "z-particle cumulant" f q (i) by 

/eawW*. (is) 



we find 



/2W = a -(b), (19) 
A(<) = a^-(b^)-2a(b)+2(b) 2 , (20) 



etc. In Section ^ below, we shall show that these f q as well as the normalization ^ orm must 
be corrected for a remaining statistical bias. This correction should become important for 
small data samples. 

Besides counting the number of combinations of q particles lying inside a sphere of ra- 
dius e, a second useful form for Star moments and cumulants are the so-called differential 
moments: Here, one defines not only a maximum distance e t but a minimum also, e t -i (t 
can define a sequence of such distances). For a given combination of q— 1 particles around 
a center particle at X^, at least one of these must lie inside the spherical shell bounded by 
radii et-i and et, while the others are restricted only by the maximum distance et. This is 
illustrated in Figure 2. 




Figure 2: Sphere counts for differentials. Given the center particle at X^, at least 
one other particle must be in the shell bounded by radii et-i and et to count. For 
q = 2, this reduces to subtracting the sphere count for et-i from that for et- Higher 
orders are also easily calculated. 



This definition leads rigorously |1[ to simple and efficient prescriptions for measurements. 
For q = 2, the unnormalized differential moment is, with A£ g (e t ) = (52^ A£ g (zi,e t )^ 

A| 2 (ii,*)=n(X? 1 ,e t )-n(X? 1 ,e t _i) = o t -a t _i, (21) 

the latter defining the shortened notation. For higher orders, we find 

Ai g (h, t) = n(Xl,e t )^ - n(Xl,et^~ 1] = ~ afc 1] , (22) 

i.e. just the difference of [q— l]th factorial powers of two sphere counts. Equivalent forms 
for the (differential) normalizations A££ orm and differential cumulants Af q are easily found, 
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leading to the normalized differential moments and cumulants 



AF q (t) 
AK q (t) 



[9-1] 



7 [9-l] 



9-1 



(23) 
(24) 



in obvious notation. The point is that these quantities can all be measured in terms of the 
two types of sphere counts, n within the same a-event and nj, within the other 6-events; see 
Eqs. PH0)- 



4 Eliminating effects of the overall multiplicity distri- 
bution 

If there are N particles within the total phase space (sky region) Q tot considered, the nor- 
malized factorial moment for this whole region is F q (tt to t) = (N^) / (N) 9 , which is unity 
only when the overall multiplicity distribution is poissonian. All measurements of F q , K q 
and their differentials thus implicitly contain correlations arising from the non-poissonian 
nature of the overall multiplicity distribution. This is as it should be, of course, but it may 
sometimes be desirable to eliminate this dependence on global effects (for example when the 
multiplicity distribution is artificial, such as in centrality cuts in heavy ion collisions). One 
way of achieving this is to modify all the formulae of the preceding sections by changing the 
event-by-event counts to 

^- k ^=jk (x - Xn)= iM^- (25) 

where N is the event's multiplicity within the total domain f2 to t> and p q to 

1 N 

h q (xi,...,x q ) = 5(xt- X h )S(x 2 - X i2 ) ■ ■ ■ 5(x q - X iq ) 

=1 

(26) 



Pq\ X li ■ ■ ■ i 



In tot p q (x 1 ,...,x q )dx l ...dx q ' 

The event average h q = (h q ) satisfies the requirements of a joint probability (normalization 
to unity and correct projection properties). These changes then propagate to yield, for 
example 

™-(i?v^) (27) 



7 



which yields F q (Q t ot) = 1 for any overall multiplicity distribution. Analogously, the individ- 
ual terms in the cumulant functions are divided by their respective integrals, so that, for 
example, 

C 3 (cci, x 2 , x 3 ) — ► cs(aji, x 2 , x 3 ) = h 3 -J2 h\h 2 + 1h x hxhi (28) 

(3) 

yielding after normalization (cf. Eq. (pC|)) 

K 3 (e) = (29) 

(^{(N a a -l)^~(jf ] )" 2 ^(^) +2 (^)))/(^?(^) )' 

Statistical independence is understood for these cumulants to mean a factorization of the 
probabilities h q rather than of the densities p q . 



5 Bose-Einstein moments and cumulants 

One great advantage of correlation integrals in general is that they allow the use of variables 
which are functions of two or more particles pi, while conventional binning is usually done 
in terms of single-particle variables only 

Bose-Einstein correlations are a prime example of the use of relative coordinates: the 
quantum mechanical interference of identical particles manifests itself in a rise of the two- 
particle correlation function 

. / „ X P2(Pl,P 2 ) i /om 

k 2{p 1 ,P 2 ) = — f \ 7 7 - 1 (30) 

Pl(Pl)Pl(P 2 ) 

at small relative momenta q = P\—p 2 - Other formulations [|7| test the correlation in terms of 
the one-dimensional variable Q 2 = —{pi—p 2 ) 2 , the relative four-momentum. The correlation 
integral formalism can be utilized for both these variables to yield moments and cumulants 
of higher order ||. Taking Q 2 as an example, one first integrates out the unneeded degrees 
of freedom in both p 2 and the normalization p\p\ 0, 

p 2 (Q 2 ) = Jd 3 p 1 d 3 p 2 p 2 ( Pl ,p 2 )5[Q 2 + (p l -p 2 ) 2 }, (31) 
Pi®Pi(Q 2 ) = Jd 3 p 1 d 3 P2 p 1 (p 1 ) Pl (p 2 )5[Q 2 + (p 1 -p 2 ) 2 ], (32) 

which, using the delta function form of Eq. translates into the measurement prescriptions 

P 2 (Q 2 ) = (E^-WM. (33) 



Am 2 _ (n ab \ 2 \ - I 

iVf 1 



Pi®Pi(g 2 ) = -p,EE^ 2 -(^) 2 ] = (E(E^ 2 -(^ 6 ) 2 ])) , (34) 

Vev aj^b i,j \ i \ j II 
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where (Qlj) 2 = —(P? — Pj) 2 measures the relative four-momentum between particles i and 
j taken from two different events a and b. Here, too, we see the direct emergence of the 
event mixing prescription as the appropriate method of normalizing Bose-Einstein correlation 
functions. 

For measurement of higher orders, one must first make an ansatz how the q three- 
momenta are to be combined into a single variable, the choice of which depends on physical 
arguments of the specific system and the signal being sought. One possibility is to sum 
all q(q—l)/2 pairs of relative four-momenta to give a measure of the overall g-particle four- 
momentum H, e.g. for q = 3, 

Q 2 = ~(pi ~ P2) 2 - (pi - P3) 2 - (P2 - p 3 ) 2 ; (35) 

this amounts to a GHP-type topology of the correlation integral in the four- momenta [l], |3|] . 
(The Q 2 defined in this way is merely the g-particle invariant mass minus a constant.) 
Moments are found by formulas analogous to Eq. (|31f) above, while cumulants are constructed 
directly from Eqs. (P^)jf. inserted into 

C q (Q 2 ) = [d 3 p 1 ...d 3 p q C q (p 11 ... lPq )5[Q 2 + J2 ( Pa -pp) 2 }, (36) 

a</3=l 

i.e. the expansion of C q in terms of the p q must be done before projection of the three- 
momenta onto Q 2 . 



6 Biased and unbiased estimators 

The use of the Star integral (or other forms such as the form used above for Bose-Einstein 
correlations) permits much more accurate measurements and hence will likely reveal more 
detailed structure of the underlying dynamics. Greater accuracy requires, however, that 
possible biases be understood on a higher level than before. One such bias arising generally 
in the measurement of correlations has to do with the theory of estimators [[UJ . 

To understand this, we must go back to the basics of sampling theory. For a given 
quantity of interest ( "statistic" ) U, there ideally exists an infinite set of measurements U ; 
this is termed the population of such measurements. A statistical average based on the whole 
population would yield the "true" value U of this quantity. 

In practice, the size of the set of measurements carried out is limited, corresponding to a 
single sample of N ev measurements taken out of the total population. Many such samples M 
could theoretically be taken, each one yielding a sample average (U), the set of which in itself 
forms a distribution, the sampling distribution. While there is no way to ascertain where 
within this distribution the (U) s obtained from a particular sample will fall, at least one can 
test whether the average of this sampling distribution coincides with U. Surprisingly, such 
a sampling average 

{U} = lim £(f/) s /A/- (37) 
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does not necessarily coincide with U except in the (for the experimentalist uninteresting) 
case Nev — > oo. When it does not, one looks for a modification, say e(U); which is called an 
unbiased estimator of U if it fulfils the condition 

{e(U)} = U for all values of N ev . (38) 

The age-old problem of finding suitable estimators has been extensively investigated and we 
merely quote the results. It has been shown that the inclusive density p q we have been using 
in the previous sections is an unbiased estimator, 

{p q (x u ...,x q )} = p q (x u ...,x q ); (39) 

in addition, we note that the sampling average of a single event inclusive density p q as defined 
in Eq. (|2]) also yields p since Eq. (|38| ) is valid for samples consisting of single events, N ev = 1, 

{p q {x 1 , ...,x q )} = p q (xx, ...,x q ). (40) 

However, whenever two or more event averages are involved, the naive product of sam- 
ple densities yields a biased estimator, for example {pi(jci)pi(a;2)} ^ p~i(x\)p~\(x2) and 
{p2(x\, x 2 )pi(x 3 )} 7^ p2(xi,x 2 )pi(x 3 ), so that all normalizations and higher-order cumu- 
lants discussed in the previous sections must be corrected to yield unbiased estimators of 
their corresponding expectation values. 

Consider for example the product of two single-particle densities. Let A be the number 
of events used for averaging; when one corrects interferometry or conventional factorial 
moments, A will be equal to N ev ; for the Star integral, the usual choice is A = N ev — 1 or, 
when a faster inner event loop is desired, A can be made much less than N ev |L]. Now, using 
Eq. ([5]), we have 

{p 1 {x l )p l {x 2 )} = {^£p1H*i)P?(* 2 )} + {^EPiH^I^EPi 2 ^)}. 

the second part being the desired true value p\{x\)p\{x2) . The reason why {pipi} does not 
yield the true value lies in the "diagonal terms" e\ = e 2 in the double sum above which 
prevent the desired factorization, as {pTpT} 7^ {fii 1 } {ft?} unless e\ and e 2 refer to two 
different (and hence independent) events. Clearly, the desired unbiased estimator is given 
by the double sum restricted to unequal events, since 

{j^EJVp?} = J[h)Z ^ = = ^ ■ < 41) 

In general, therefore, the unbiased estimator for the product of q single-particle densities is 
given by 

^ £ (42) 

ei^e 2 ^...^eq 
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From this prescription, we obtain the unbiased estimators for the Star integral normalization 
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Pnorm 
S2 


= (b), 






Pnorm 
S3 


= (b) 2 ~ 


K 2 (b, b) 

{A-iy 




Pnorm 
S4 


= (bf- 


3(b)K 2 (b,b) 

(A-l) 


2K 3 (b,b,b) 

(A - 1)P] ' 


Pnorm 
S5 


= (b) 4 ~ 


Q(b) 2 K 2 (b,b) 

(A-l) 


8{b)K 3 {b, b, b) + 3k 2 (&, b) 
(A-l)! 2 ! 



(43) 
(44) 

(45) 



6K 4 (b,b,b } b) +9K 2 2 (b, b) 
(A-l) PI 



(46) 



where 



k 2 (U,V) = (UV)-(U)(V), (47) 
k 3 (U,V,W) = (UVW)-J2(UV)(W)+2(U)(V)(W), (48) 

(3) 

Ki(u,v,w,x) = (uvwx) -J2(u)(vwx) -J2(uv)(wx) 

(4) (3) 

+ 2J2(U)(V)(WX) - 6(U)(V)(W)(X) , (49) 

(6) 

where the sums indicate the number of combinations to be taken and U, V, ... is any statistic 
of interest; for example n 2 (b,b^) = (bb^) — (b)(b^). Note that (for the Star integral) the 
second order normalization does not need a correction; this is because the first event sum 
over ei is always pulled out in calculating the a-quantities; only the sums over e 2 , e 3 . . . must 
be made explicitly unequal. 

For the cumulants, we need the more general statement: if p qi , p q2 . . . p qK are densities 
of order q 1 , q 2 , . . . q K , the unbiased estimator of their product is given by 



AW 

Implementing these, we find the unbiased estimators for the z-particle cumulants to be 

hit) = a -(b), (51) 
/»(<) = aW-(bM)-2a(b)+2(b) 2 - J ?- 1 K 2 (b,b), (52) 
/ 4 (i) = a [3] -(b [3] }-3a [2] (b}-3a(b [2] } + Q(b}(b [2] }+Qa(b} 2 -Q(b} 3 

+ A~^l i {3{b) ~ a) ^ h) ~ ^ bl2]) } ~ (A-l)W K * {b ' 6 ' h) ' (53) 
/ 5 (i) = a® - (6 [4] ) -4aW(b) - Aa(b [3] ) 
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- 6a® (b®) + 8(b) (b®) + 12^ (6> 2 + Q(b^}(b^} 
+ 24a(6)(6 [2] ) - 36(6) 2 (6 [2] ) - 24a(6) 3 + 24(6) 4 

- -r^— {(6a [2] - 18(6 [2] ) - 36a(6) + 72(6) 2 ) k 2 (6, b) 

J\ — 1 

+4k 2 (6,6^) +3K 2 (bM,bM) + (12a- 36(6)) k 2 (6,6 [2] )} 

24 

+ (A _ 1)[2] 6) + (8(6) - 2a) K3 (6, 6, b) - 3K 3 (b, 6, b^)} 

72 

- {A _ 1)[3] {2^4(6,6,6,6) + 3^(6,6)} . (54) 

For very small samples, when the inner event average sum J2b cannot be taken strictly 
over 6 7^ a, corrections must also be made for equal-event terms JT(J. These are very 



important for small samples found e.g. in fixed- iV cuts and cosmic ray data. When full event 
mixing is implemented for Bose-Einstein correlations and conventional factorial moments or 
cumulants, similar bias corrections are mandatory ||10||. 
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